18.1 Crime and punishment

The data in linearRegression_crimePunishment.csv contains the murder rate per capita and the rate of automobile crimes per 100,000 individuals (both on the log scale) in the ten US States that have changed their legislation on capital punishment since 1960 (in all cases the states abolished capital punishment). We also include a dummy variable (“law”) that is 1 if the state allows capital punishment in that year, and 0 otherwise. The crime data is from http://www.disastercenter.com.

Problem 18.1.1. Graph the data and comment on any trends.

crime_data <- read.csv('linearRegression_crimePunishment.csv')
crime_data

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ────────────────────────────────────────────────────── tidyverse 1.3.0 ──
✓ ggplot2 3.3.2     ✓ purrr   0.3.4
✓ tibble  3.0.3     ✓ dplyr   1.0.1
✓ tidyr   1.1.1     ✓ stringr 1.4.0
✓ readr   1.3.1     ✓ forcats 0.5.0
── Conflicts ───────────────────────────────────────────────────────── tidyverse_conflicts() ──
x purrr::accumulate() masks foreach::accumulate()
x dplyr::filter()     masks stats::filter()
x dplyr::lag()        masks stats::lag()
x purrr::when()       masks foreach::when()

ggplot(data=crime_data) + geom_point(mapping=aes(x=year, y=murder))

ggplot(data=crime_data) + geom_line(mapping=aes(x=year, y=murder))

ggplot(data=crime_data) + geom_smooth(mapping=aes(x=year, y=murder,linetype=stateName))

ggplot(data=crime_data) + geom_smooth(mapping=aes(x=year, y=murder,group=stateName))

ggplot(data=crime_data) + geom_smooth(mapping=aes(x=year, y=murder, color=stateName))

#ggplot(data=filter(crime_data, stateName=='Connecticut')) + geom_line(mapping=aes(x=year, y=murder)) + geom_line(mapping=aes(x=year, y=car))

ggplot(data=filter(crime_data, stateName=='Connecticut'), aes(year)) + geom_line(aes(y=murder, color="murder")) + geom_line(aes(y=car/2, color="car")) + geom_line(aes(y=law, color='law'))

filter(crime_data, stateName=='Connecticut')

df <- filter(crime_data, stateName=='Connecticut') %>%

Warning messages:
1: In readChar(file, size, TRUE) : truncating string with embedded nuls
2: In readChar(file, size, TRUE) : truncating string with embedded nuls
3: In readChar(file, size, TRUE) : truncating string with embedded nuls
4: In readChar(file, size, TRUE) : truncating string with embedded nuls
5: In readChar(file, size, TRUE) : truncating string with embedded nuls
6: In readChar(file, size, TRUE) : truncating string with embedded nuls
7: In readChar(file, size, TRUE) : truncating string with embedded nuls
8: In readChar(file, size, TRUE) : truncating string with embedded nuls

  select(year, murder, car, law) %>%
  gather(key="variable", value = "value", -year)

ggplot(df, aes(x=year, y=value)) +
  geom_line(aes(color=variable, linetype=variable)) + ggtitle("Connecticut") + theme(plot.title = element_text(hjust = 0.5)) 

NA
NA

Warning messages:
1: In readChar(file, size, TRUE) : truncating string with embedded nuls
2: In readChar(file, size, TRUE) : truncating string with embedded nuls
3: In readChar(file, size, TRUE) : truncating string with embedded nuls
4: In readChar(file, size, TRUE) : truncating string with embedded nuls
5: In readChar(file, size, TRUE) : truncating string with embedded nuls
6: In readChar(file, size, TRUE) : truncating string with embedded nuls
7: In readChar(file, size, TRUE) : truncating string with embedded nuls

for (lstate in unique(crime_data$stateName)){
  df <- filter(crime_data, crime_data$stateName == lstate) %>%
  select(year, murder, car, law) %>%
  gather(key="variable", value = "value", -year)
  
  print(ggplot(df, aes(x=year, y=value)) +
  geom_line(aes(color=variable, linetype=variable)) +  ggtitle(lstate) + theme(plot.title = element_text(hjust = 0.5)))
  
}

Problem 18.1.2.

A simple model for murder rates is of the form,

        murderi,t ∼ N (α + βpenaltyi,t + γcari,t, σ) (18.1)

where we assume that the effect of having the death penalty is given by β, which is assumed to be the same across all states. We include cari,t – a measure of crimes on automobiles, as an independent variable to proxy for the contemporaneous underlying level of crime. Estimate this model and hence determine whether the death penalty acts as a deterrent to murder.

library(rstan)
options(mc.cores = parallel::detectCores())

dataList <- list(N = 540, murder = crime_data$murder, car = crime_data$car,

There were 16 warnings (use warnings() to see them)

                 law = crime_data$law, state = crime_data$state)

fit <- stan(file = 'Model_Problem_18_1.stan', data = dataList, iter=1000, chains=4, seed=42)

starting worker pid=58687 on localhost:11242 at 17:39:39.666
starting worker pid=58701 on localhost:11242 at 17:39:39.973
starting worker pid=58715 on localhost:11242 at 17:39:40.713
starting worker pid=58734 on localhost:11242 at 17:39:41.209

SAMPLING FOR MODEL 'Model_Problem_18_1' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000133 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.33 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'Model_Problem_18_1' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 0.000158 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.58 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:   1 / 1000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'Model_Problem_18_1' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 0.000174 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.74 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)

SAMPLING FOR MODEL 'Model_Problem_18_1' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 0.000169 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.69 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 3: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 2: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 4: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 3: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 2: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 4: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 3: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 2: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 4: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 3: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 2: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 4: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 3: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 3: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 2: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 2: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 4: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 4: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 3: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 2: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 4: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 3: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 4: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 2: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 3: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 4: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 2: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 3: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 1.16125 seconds (Warm-up)
Chain 1:                1.0339 seconds (Sampling)
Chain 1:                2.19514 seconds (Total)
Chain 1: 
Chain 4: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 3: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 1.13128 seconds (Warm-up)
Chain 3:                1.04703 seconds (Sampling)
Chain 3:                2.17831 seconds (Total)
Chain 3: 
Chain 2: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 4: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 1.08466 seconds (Warm-up)
Chain 4:                1.16152 seconds (Sampling)
Chain 4:                2.24619 seconds (Total)
Chain 4: 
Chain 2: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 1.27196 seconds (Warm-up)
Chain 2:                1.17677 seconds (Sampling)
Chain 2:                2.44873 seconds (Total)
Chain 2: 

Warning messages:
1: In readChar(file, size, TRUE) : truncating string with embedded nuls
2: In readChar(file, size, TRUE) : truncating string with embedded nuls
3: In readChar(file, size, TRUE) : truncating string with embedded nuls
4: In readChar(file, size, TRUE) : truncating string with embedded nuls
5: In readChar(file, size, TRUE) : truncating string with embedded nuls
6: In readChar(file, size, TRUE) : truncating string with embedded nuls
7: In readChar(file, size, TRUE) : truncating string with embedded nuls
8: In readChar(file, size, TRUE) : truncating string with embedded nuls

print(fit)

Inference for Stan model: Model_Problem_18_1.
4 chains, each with iter=1000; warmup=500; thin=1; 
post-warmup draws per chain=500, total post-warmup draws=2000.

        mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
alpha  -1.48    0.01 0.22  -1.91  -1.63  -1.48  -1.34  -1.05   784    1
beta    0.24    0.00 0.06   0.12   0.20   0.24   0.28   0.35  1170    1
gamma   0.45    0.00 0.04   0.37   0.43   0.45   0.48   0.53   781    1
sigma   0.62    0.00 0.02   0.59   0.61   0.62   0.64   0.66  1200    1
lp__  -16.05    0.05 1.32 -19.24 -16.77 -15.77 -15.06 -14.38   826    1

Samples were drawn using NUTS(diag_e) at Tue Sep  8 17:39:58 2020.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).